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Fourier  spectral  method  can  achieve  exponential  accuracy  both  on  the  approximation  level 
and  for  solving  partial  differential  equations  if  the  solutions  are  analytic.  For  a  linear  partial 
differential  equation  with  a  discontinuous  solution,  Fourier  spectral  method  produces  poor  point- 
wise  accuracy  without  post-processing,  but  still  maintains  exponential  accuracy  for  all  moments 
against  analytic  functions.  In  this  note  we  assess  the  accuracy  of  Fourier  spectral  method  applied 
to  nonlinear  conservation  laws  through  a  numerical  case  study.  We  find  that  the  moments  with 
respect  to  analytic  functions  are  no  longer  very  accurate.  However  the  numerical  solution  does 
contain  accurate  information  which  can  be  extracted  by  a  post-processing  based  on  Gegenbauer 
polynomials. 
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1  Introduction 


In  this  note  we  are  concerned  with  the  accuracy  of  Fourier  spectral  method  for  the  solution  of  a 
nonlinear  conservation  law 

dtu  +  dxf(u)  -  0,  - 1  <  x  <  1 

u(z,0)  =  u°(x)  (1.1) 

where  the  initial  condition  u°(x)  is  2-periodic.  As  is  well  known,  solutions  to  (1.1)  typically  con¬ 
tain  discontinuities  even  if  the  initial  condition  u°(x)  is  analytic  (In  this  paper,  for  simplicity  of 
presentations,  we  will  use  analytic  functions  to  represent  general  smooth  functions;  similar  results 
can  also  be  obtained  for  Ck  or  C°°  functions).  The  purpose  of  this  note  is  to  assess  accuracy  under 
such  a  situation  through  a  numerical  study. 

We  start  by  recalling  the  Fourier  approximation  operator  Sn  to  an  L2  function  w(x): 

N 

Snu(x)  =  £  ukeik™  (1.2) 

ks-N 

where  the  Fourier  coefficients  uk  are  defined  by 


uk  =  \j\*{x)e-ik™dx 


for  Fourier  Galerkin,  and  by 


1  v  ,  .  —ikirx  2,7 

=  2777  £  “(’i)e  *’  *’  ~  2777 


J=-N 


for  Fourier  collocation.  We  will  also  use  the  notation  uh(x)  =  S^u(x).  To  solve  the  partial 
differential  equation  (1.1),  the  standard  Fourier  spectral  algorithm  is 

SN(dtVN  +  dxf(vN))  =  0,  -1  <  x  <  1 

u^(x,0)  =  Stvu°(x)  (1-5) 


where  vh(z,  t)  =  J2k=-N  Vk(t)e'k*x  is  supposed  to  approximate  the  exact  solution  u(x,t)  of  (1.1), 
and  Sn  is  either  the  Galerkin  or  the  collocation  Fourier  approximation  operator  defined  by  (1.2)- 
(1.3)  or  by  (1.2)-(1.4). 


The  approximation  error 


tt(i)  -  Snu(x)  (1.6) 

is  well  known  to  be  exponentially  small  (i.e.,  it  is  O(r^)  for  0  <  r  <  1)  if  u(x)  is  analytic.  However, 
if  u(x)  is  only  piecewise  analytic  and  discontinuous,  the  approximation  error  (1.6)  is  0(1)  near  the 
discontinuity  and  O(^)  elsewhere.  This  is  known  as  the  Gibbs  phenomenon  (see,  e.g.,  [4]  and  [3]). 
Fortunately,  even  if  the  accuracy  is  poor  in  the  point-wise  sense,  it  is  still  excellent  for  the  moments 
against  any  analytic  functions.  For  any  L 2  function  u(x)  and  any  analytic  function  w(x),  we  have 
[5]: 


J  (u (x)  -  u\(xs'\w(x)dx 


<  Cr 


N 


(1.7) 


for  some  constant  C  and  0  <  r  <  1.  This  property  is  the  \osis  of  all  the  “reconstruction”  or  '  post¬ 
processing”  techniques.  These  techniques  try  to  recover  expcnen'ial  or  at  least  high  order  accuracy 
for  point  values  based  on  the  Fourier  approximation  Snu(x)  of  a  piecewise  analytic  function.  In 
other  words,  one  tries  to  obtain  a  small  post-processed  approximation  error 


«(*)  ~  PnSnu(x) 


(1.8) 


where  Ppj  is  some  post-processing  operator.  Examples  of  Pjy  include  various  high  frequency  filters 
[14],  [11],  [16],  [2],  which  are  of  the  form 

N  •  /fc\ 

PNSNu(*)=  £  (1.9) 


with  Snu(x)  given  by  (1.2).  The  function  <t(£)  in  (1.9)  is  even  (or  satisfies  <r(-£)  =  <r(f)  if  it 
is  complex  valued  as  in  [2]),  smooth  (the  accuracy  of  the  filter  depends  upon  its  smoothness), 
supported  in  (—1, 1)  and  satisfies  cr(O)  =  1  and  a^k\0)  =  0  for  1  <  k  <  K  (with  accuracy  of  the 
filter  again  depending  on  K).  These  filters  can  recover  high  order  or  even  exponential  accuracy 
in  the  smooth  regions  away  from  the  discontinuities  (the  filter  in  [2]  can  also  recover  high  order 
accuracy  up  to  the  discontinuity  from  one  side).  A  more  recent  example  of  Pn  is  the  Gegenbauer 
polynomial  based  procedure  discussed  in  [6],  [7],  [8],  [9]  and  [10],  which  can  give  uniform  exponential 
accuracy  for  all  x  right  up  to  the  discontinuity  for  piecewise  analytic  functions.  In  this  sense  spectral 
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Fourier  approximation  is  also  exponentially  accurate  for  piecewise  analytic  functions — one  only  has 
to  extract  the  hidden  information  from  the  poor  approximation  S/\r(x)  using  the  post-processor 
Pn- 

When  spectral  method  is  used  to  solve  the  PDE  (1.1),  we  can  consider  the  following  different 
types  of  errors.  The  strongest  is  the  point- wise  error  from  the  exact  solution  u(x,t): 

u(x,t)~  vN(x,t),  (1.10) 

which  cannot  be  small  even  for  t  =  0  due  to  the  Gibbs  phenomenon.  A  more  reasonable  error  is 

the  point-wise  error  of  the  numerical  solution  «jv(x,t)  from  the  Fourier  approximation  of  the  exact 
solution  u/v(z,t): 

UN(x,t)  -  vN(x,t).  (1-11) 

If  this  error  is  exponentially  small,  we  claim  the  spectral  method  for  (1.1)  is  exponentially  accurate 
because  of  the  post-processor  (1.8)  for  the  exact  solution  u(x,t).  An  even  weaker  error  is  defined 
by  the  error  in  the  first  few  Fourier  coefficients,  i.e. 

Uk(t)-Vk(t)  (1.12) 

for  the  first  few  k,  where  Uk{t)  are  the  Fourier  coefficients  of  the  exact  solution  u(x,t)  of  (1.1). 
This  is  actually  an  example  of  the  more  general  definition  of  error  in  moments  with  respect  to  any 
analytic  function  w(x): 

J  (un(x)  —  vn(x))w(x)<1x  (1.13) 

In  fact,  as  long  as  this  error  in  moments  is  exponentially  small,  we  claim  that  the  spectral  method 
is  exponentially  accurate  in  solving  (1.1)  by  using  property  (1.7)  for  the  exact  solution  u(a:,t)  and 
the  post-processing  (1.8). 

If  the  PDE  (1.1)  is  linear  (i.e.  f(u)  =  a(x,t)u ),  it  is  proven  in  [5],  [1]  that  spectral  Fourier 
method  is  exponentially  accurate  in  the  sense  that  (1.13)  is  exponentially  small.  A  post-processing 
(1.8)  applied  to  vn(x,  t )  would  then  yield  an  exponentially  accurate  point-wise  approximation  to  the 
exact  solution  u(x,t).  However,  if  (1.1)  is  nonlinear,  it  is  still  a  theoretically  open  problem  whether 
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spectral  Fourier  method,  equipped  with  either  high  frequency  filtering  or  vanishing  viscosity  [15], 

[12] ,  is  exponentially  (or  high  order)  accurate  in  the  sense  of  (1.13).  Computational  evidence  in 

[13]  seems  to  suggest  that,  even  in  this  nonlinear  case,  highly  accurate  information  is  still  implicitly 
contained  in  the  numerical  solution  and  can  be  extracted  (at  least  away  from  the  discontinuity) 
by  a  post-processing  using  high  frequency  filtering.  In  the  next  section  we  will  perform  a  detailed 
numerical  case  study  about  this  accuracy  issue  for  Burgers’  equation  (/(u)  =  ^-).  We  use  a 
high  frequency  solution  filter  to  stablize  the  algorithm,  and  post  process  the  numerical  result 
using  the  procedure  based  Gegenbauer  polynomials  [6],  [10].  We  observe  that  the  spectral  Fourier 
method  is  not  very  accurate  in  the  sense  of  moments  against  analytic  functions  (1.13).  However, 
numerical  evidence  does  indicate  the  possibility  of  very  high  accuracy  under  some  weaker  definition 
of  accuracy,  perhaps  some  average  of  Fourier  coefficients,  since  the  post-processed  result  Pnvm(x,  t) 
is  much  more  accurate  than  the  Fourier  coefficients  themselves,  and  accurate  Fourier  coefficients 
can  be  “reconstructed”  from  this  post-processed  solution  Pnvn(x,  t). 

2  A  Numerical  Case  Study  about  Accuracy 

In  the  numerical  solution  reported  in  this  section,  time  discretization  is  by  a  third  order  Runge- 
Kutta  method,  with  a  time  step  At  sufficiently  small  such  that  the  spatial  error  is  dominant 
in  all  cases.  We  compute  the  exact  solutions  of  the  PDE  by  Newton  iterations  on  the  implicit 
characteristic  equations,  and  compute  the  Fourier  coefficients  of  a  function  (if  not  analytically 
given)  by  using  a  sufficiently  accurate  numerical  quadrature. 

We  first  solve  a  linear  equation 

3 

QtU  +  - - - - -  QxU  —  0,  -1  <  X  <  1 

5-4cos(7rx) 

u(x,0)  =  x  (2.1) 

with  periodic  boundary  conditions,  up  to  t=l,  using  the  Fourier  Galerkin  method: 

Sn  (dtvw  +  - — — -  - — -  dxvN  )  =  0, 

V  5-4cos(?rx)  / 
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vN(x,0) 


Snx  = 


(2.2) 


k  =  -N 
k  ^  0 

Standard  Galerkin  method  is  stable  for  this  linear  problem  but  produces  poor  point  value  accuracy 
(Figure  1,  left).  However,  the  accuracy  in  the  first  few  Fourier  coefficients,  as  representatives  of 
moments  against  analytic  functions,  are  computed  more  accurately  (Figure  1,  right). 


Fig.  1:  Errors  in  log  scale,  linear  PDE  (2.1).  Fourier  Galerkin  using  2 N  +  1  modes,  for  N  =  10; 
N  =  20;  N  —  40  and  N  =  80.  Left:  point-wise  errors;  Right:  errors  in  the  first  10  Fourier 
coefficients. 

In  order  to  compare  with  the  nonlinear  case  reported  later,  we  solve  the  same  linear  equation 
(2.1)  using  the  filtered  Fourier  method,  i.e.,  after  each  Runge-Kutta  time  step,  the  numerical 
solution  is  filtered  by  (1.9)  with  the  exponential  filter: 

<x(0  =  e-“Klr  (2.3) 

where  r  is  increasing  with  N  and  is  related  to  the  order  of  the  filter,  and  a  is  chosen  such  that  e~a 
equals  machine  zero  (  10-16  for  double  precision).  The  exponential  filter  (2.3)  has  the  advantage 
of  simplicity,  and  usually  it  works  equally  well  as  more  complicated  filters  [16].  For  this  linear 
problem,  as  well  as  for  the  nonlinear  Burgers’  equation  later,  we  will  use  the  Fourier  method  with 
the  following  choice  of  filter  orders:  r  =  4  for  IV  =  10;  r  =  6  for  N  =  20;  r  =  8  for  AT  =  40  and 
r  =  12  for  N  =  80.  The  result  is  shown  in  Figure  2.  Comparing  with  Figure  1,  we  can  see  better 
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point  value  accuracy  in  the  smooth  region  because  of  the  filters,  and  similar  (good)  accuracy  for 
the  first  few  Fourier  coefficients. 


Fig.  2:  Errors  in  log  scale,  linear  PDE  (2.1).  Fourier  Galerkin  using  27V  +  1  modes  with  exponential 
solution  filters  of  order  r.  r  =  4  for  TV  =  10;  r  =  6  for  TV  =  20;  r  =  8  for  TV  =  40  and  r  =  12  for 
TV  =  80.  Left:  point-wise  errors;  Right:  errors  in  the  first  10  Fourier  coefficients. 

The  computational  result  for  the  linear  equation  is  not  surprising  since  it  just  shows  the  proven 
fact  [5],  [1]  that  Fourier  coefficients,  as  representatives  of  moments  against  analytic  functions,  are 
computed  with  exponential  accuracy  by  the  spectral  Fourier  method,  and  filtering  will  recover 
exponential  point  value  accuracy  in  smooth  regions  away  from  the  discontinuity.  It  should  be 
noticed  that,  for  the  same  TV,  the  accuracy  for  the  first  few  Fourier  coefficients  is  at  the  same  level 
or  better  than  the  best  point  value  accuracy  in  the  smooth  region  after  filtering.  This  is  again  not 
surprising  since  point  value  accuracy  is  obtained  from  the  Fourier  coefficients  through  filtering. 

We  now  consider  the  nonlinear  problem  we  are  really  interested  in:  we  solve  the  nonlinear 
Burgers’  equation 

dtu  +  dx  ^ y 

u(x,0)  =  0.3  +  0.7  sin(irx).  (2.4) 

The  solution  develops  a  shock  at  t  =  and  we  compute  the  solution  up  to  t  =  1.  The  initial 
condition  is  chosen  such  that  the  shock  is  moving  with  time.  For  this  nonlinear  PDE,  the  standard 


j  =0,  -1  <  x  <  1 
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Galerkln  method  cannot  converge  to  the  entropy  solution  [15].  One  would  need  to  add  dissipations 
either  by  the  high  frequency  solution  filtering  (1.9)  or  by  the  spectral  vanishing  viscosity  [15],  [12], 
[13].  Numerical  results  for  the  Burgers’  equation  with  the  vanishing  viscosity  method  can  be  found 
in,  e.g.,  [13].  Here  we  will  only  report  the  results  obtained  by  solution  filtering,  using  the  same  r 
as  in  the  previous  linear  case  (2.1).  We  have  also  computed  with  the  vanishing  viscosity  methods 
and  have  obtained  similar  results. 

In  Figure  3  we  plot  the  point-wise  error  u(z,t)  —  vs{x,t)  (left),  and  the  error  for  the  first  10 
Fourier  coefficients  (right).  While  the  pattern  of  the  point-wise  errors  are  similar  to  the  linear  case 
in  Figure  2,  the  errors  for  the  Fourier  coefficients  are  clearly  much  worse  in  comparison.  As  a  matter 
of  fact,  for  the  same  N,  the  errors  for  the  first  few  Fourier  coefficients  are  a  few  magnitudes  larger 
than  the  smallest  point  value  error  in  the  smooth  region.  This  is  clearly  different  from  what  we 
observe  in  the  linear  case  in  Figure  2,  and  suggests  that  the  first  few  Fourier  coefficients,  again  as 
representatives  of  moments  against  analytical  functions,  are  no  longer  computed  with  exponential  or 
high  order  accuracy.  It  is  sort  of  puzzling  that  each  difference  in  the  Fourier  coefficients  «*(<)- v*(t) 
is  relatively  large  (Figure  3  right),  but  the  point-wise  error  u(x,t)-~VN(x,t),  which  is  just  an  average 
of  u;fc(t)  -  ik{t)  (against  ar.  0(1)  function  etknx),  is  much  smaller  in  the  smooth  region  (Figure  3 
left).  Clearly  some  cancellation  is  present. 


Fig.  3:  Errors  in  log  scale,  Burgers  equation  (2.4).  Fourier  Galerkin  using  2 N  -I-  1  modes  with 
exponential  solution  filters  of  order  r.  r  =  4  for  N  =  10;  r  =  6  for  N  =  20;  r  =  8  for  N  =  40  and 
r  =  12  for  N  =  80.  Left:  point-wise  errors;  Right:  errors  in  the  first  10  Fourier  coefficients. 
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Next,  we  apply  the  Gegenbauer  post-processor  [6]  to  vs{x,t).  This  procedure  can  be  roughly 
described  as  follows:  given  the  Fourier  partial  sum  uyv(x)  of  an  analytic  but  not  periodic  function 
u(x),  one  first  finds  the  approximations  to  the  first  m  Gegenbauer  expansion  coefficients  of  the 
function  u(x).  Here  Gegenbauer  polynomials  are  orthogonal  polynomials  in  [-1, 1]  under  the  weight 
function  (1  —  i2)A-a .  One  then  uses  this  Gegenbauer  series  with  the  first  m  terms  to  approximate 
u(x)  everywhere  in  [-1, 1].  To  use  this  procedure,  one  must  know  the  location  of  the  discontinuity 
(however,  the  procedures  in  [8]  allows  one  to  handle  the  case  where  the  location  is  not  known 
exactly),  and  to  choose  the  parameters  X  and  m.  It  is  proven  in  [6]  that  when  rn  and  X  are  both 
chosen  proportional  to  (but  less  than)  JV,  the  reconstructed  point  values  are  exponentially  accurate 
everywhere  inside  [-1,1].  Thus  Gibbs  phenomenon  is  completely  removed.  The  details  can  be 
found  in  [6],  [7],  [8],  [9]  and  [10]. 

We  would  like  to  point  out  that  there  is  no  theoretical  justification  in  doing  this  post  processing 
for  the  current  nonlinear  case,  since  the  post-processing  procedure  assumes  that  the  Fourier  coef¬ 
ficients  are  accurate,  which  is  not  true  any  more.  However,  the  post-processed  result  is  surpris¬ 
ingly  good  (Figure  4).  We  can  observe  good  accuracy  everywhere  including  at  the  discontinuity 
x  —  ±1  +  0.3  as  in  the  approximation  test  case  discussed  [6].  The  reconstructed  Fourier  coef¬ 
ficients,  namely  the  Fourier  coefficients  of  P^v^f(x,  t),  are  much  more  accurate  than  before  the 
post-processing  (compare  Figure  4  right  with  Figure  3  right). 


Fig.  4:  Errors  in  log  scale,  Burgers  equation  (2.4).  Fourier  Galerkin  using  2 N  -f  1  modes  with 
exponential  solution  filters  of  order  r.  r  =  4  for  IV  =  10;  r  =  6  for  N  =  20;  r  =  8  for  IV  =  40 
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and  r  =  12  for  N  =  80.  Gegenbauer  post -processed,  with  parameters  A  =  2,m  =  1  for  N  =  10; 
A  =  3,m  =  3  for  N  =  20;  A  =  12,  m  =  7  for  N  =  40  and  A  =  62,  m  —  15,  for  N  —  80.  Left: 
point-wise  errors;  Right:  errors  in  the  first  10  Fourier  coefficients. 

This  suggests  that,  '■ven  if  viv(x,t)  or  its  Fourier  coefficients  i>k(t)  are  not  very  accurate,  it 
contains  accurate  information  which  is  extracted  in  this  case  by  the  Gegenbauer  polynomial  based 
post-processor  P#.  This  numerical  evidence  suggests  that  in  the  nonlinear  PDE  case,  Fourier 
coefficients  ik{t),  just  like  point-wise  values  in  the  linear  (or  nonlinear)  PDE  case,  are  no  longer  good 
indicators  of  accuracy.  They  themselves  are  not  very  accurate,  but  they  implicitly  contain  accurate 
information  which  can  be  extracted  by  adequate  post-processors  P/v.  This  accurate  information 
might  be  contained  in  some  averages  of  the  Fourier  coefficients  (since  the  post- processing  procedure 
based  on  Gegenbauer  polynomials  [6]  uses  certain  averages  of  Fourier  coefficients  rather  than  the 
coefficients  themselves). 

We  finally  make  two  remarks: 

Remark  2.1.  In  the  Gegenbauer  reconstruction  procedure  above  we  have  used  the  exact  shock 
location.  The  procedure  in  [8]  allows  us  to  use  an  approximate  shock  location,  determined  from 
the  Fourier  coefficients  themselves  (e.g.,  [2]).  Similarly  good  results  can  be  obtained  when  the 
reconstruction  is  performed  in  a  slightly  smaller  sub-interval  inside  which  the  solution  is  guaranteed 
to  be  analytic.  For  example,  we  use  the  shock  detector  in  [2],  which  in  this  case  detects  the  shock 
location  to  within  0.0000025  for  all  the  N  values  used,  and  a  reconstruction  inside  the  sub-interval 
[-0.999997,0.999997],  which  is  just  slightly  smaller  than  [-1, 1]  (when  numerically  detected  shock 
is  shifted  to  a;  =  —1)  and  guarantees  that  the  true  shock  is  outside  this  region.  The  result  is  shown 
in  Figure  5  (left).  It  is  clearly  as  good  as  the  case  where  the  exact  shock  location  is  used  (compare 
with  Figure  4  left). 

Remark  2.2.  If  we  use  collocation  (1.4)  instead  of  Galerkin,  (for  the  reconstruction  procedure, 
see  [10]),  the  result  is  almost  identically  good:  Compare  Figure  5  (right)  with  Figure  4  (left). 
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Fig.  5:  Point-wise  Errors  in  log  scale,  Burgers  equation  (2.4).  Fourier  method  using  2 N  +  1  modes 
with  exponential  solution  filters  of  order  r.  r  =  4  for  N  =  10;  r  =  6  for  N  =  20;  r  =  8  for  N  =  40 
and  r  =  12  for  N  =  80.  Lf  Galerkin,  Gegenbauer  post-processed  with  a  numerically  determined 
shock  location  using  the  techniques  in  [2],  which  for  this  problem  produce  shock  locations  to  within 
0.0000025  for  all  the  N  used.  The  reconstruction  sub-interval  is  [-0.999997,0.999997]  when  the 
numerical  shock  is  shifte  *  to  x  =  -1.  Parameters:  A  =  2,  m  =  1  for  N  =  10;  A  =  3,m  =  3 
for  N  =  20;  A  =  26,  m  =  9  for  /V  =  40  and  A  =  52,  m  =  17,  for  N  =  80.  Right:  collocation. 
Gegenbauer  post-processed,  with  parameters  A  =  2,  m  =  1  for  N  =  10;  A  =  3,  m  -  3  for  N  =  20; 
A  =  26,  rn  =  9  for  N  -  40  and  A  =  60,  m  =  15,  for  N  =  80. 

3  Concluding  Remarks 

Through  a  careful  numerical  case  study  for  the  Burgers’  equation,  we  have  found  that  the  Fourier 
spectral  method,  equipped  with  spectrally  small  dissipations  in  the  form  of  high  frequency  filters 
or  vanishing  viscosities,  are  not  accurate  in  the  first  few  Fourier  coefficients,  or  in  moments  against 
smooth  functions.  However,  accurate  information  is  indeed  contained  in  the  numerical  solution, 
and  can  be  extracted  by  using  the  Gegenbauer  polynomial  based  post-processor  in  [6]-[10]. 
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